Preliminaries

# load packages
library(tidyverse)
library(ggplot2)
library(sf) # simple features - store & manipulate spatial data
library(reshape2)
library(gridExtra)

# set seed
set.seed(2024)

Setup Intervention & Spillover Polygons

### Create intervention & spillover polygons ###
# Define the intervention zone polygon (333,333) to (666,666)
intervention_zone <- st_polygon(list(matrix(c(333, 333,
                                              666, 333,
                                              666, 666,
                                              333, 666,
                                              333, 333), 
                                            ncol = 2, byrow = TRUE)))

# Define the spillover zone polygon (133,133) to (866,866)
spillover_zone <- st_polygon(list(matrix(c(133, 133,
                                           866, 133,
                                           866, 866,
                                           133, 866,
                                           133, 133), 
                                         ncol = 2, byrow = TRUE)))

# Create sf objects for both polygons (simple feature collection of geometries) (CRS=NA is unit-based instead of map-based)
intervention_sf <- st_sfc(intervention_zone, crs = NA_crs_)
spillover_sf <- st_sfc(spillover_zone, crs = NA_crs_)
### Plot polygons to check ###
# Create simple data frames for both polygons
intervention_zone_df <- st_as_sf(data.frame(geometry = intervention_sf, zone = "Intervention"))
spillover_zone_df <- st_as_sf(data.frame(geometry = spillover_sf, zone = "Spillover"))

# Combine both data frames into one for plotting
zones_df <- rbind(intervention_zone_df, spillover_zone_df)

# Plot using ggplot2
ggplot() +
  # Plot spillover zone in blue
  geom_sf(data = spillover_zone_df, fill = "blue", alpha = 0.3, color = "black", size = 1) +
  # Plot intervention zone in red
  geom_sf(data = intervention_zone_df, fill = "red", alpha = 0.6, color = "black", size = 1) +
  # Set limits to match the grid size
  coord_sf(xlim = c(0, 1000), ylim = c(0, 1000), expand = FALSE) +
  # set plot theme
  theme_minimal() +
  # assign title/x/y labels
  labs(title = "Intervention and Spillover Zones w/o simulated data", x = "X Coordinate", y = "Y Coordinate")

Simulate Data

### Simulate coordinates for Cells ###
# Create a 1000x1000 grid
grid_size <- 1000

# Simulate 2000 random points on the grid with integer coordinates
points <- tibble(x = sample(1:grid_size, 2000, replace = TRUE), 
                 y = sample(1:grid_size, 2000, replace = TRUE))

# Convert points to sf (simple features) objects
points_sf <- st_as_sf(points, coords = c("x", "y"), crs = NA_crs_)

# Compute the distance from each point to the intervention zone & add to points df
distances <- st_distance(points_sf, intervention_sf)
points$intervention_distance <- as.vector(distances)

# Add indicators for if point is inside the intervention zone (1 if inside, 0 if not)
points$intervention_zone <- ifelse(st_within(points_sf, intervention_sf, sparse = FALSE), 1, 0)

# Add indicators for if point is inside the spillover zone (outside intervention zone)
points$spillover_zone <- ifelse(
  st_within(points_sf, spillover_sf, sparse = FALSE) & !st_within(points_sf, intervention_sf, sparse = FALSE), 1, 0)

# Add column "zone" as intervention/spillover/outside for each point
points$zone <- ifelse(points$intervention_zone == 1, "Intervention",
                      ifelse(points$spillover_zone == 1, "Spillover", "Outside"))
### Plot points on top of grid colored by zone ###
# Convert the points into an sf object using the coordinates (x, y)
points_sf <- st_as_sf(points, coords = c("x", "y"), crs = NA_crs_)

# Plot the polygons and the points
ggplot() +
  # Plot the spillover zone in blue
  geom_sf(data = spillover_sf, fill = "blue", alpha = 0.2, color = "black", size = 1) +
  # Plot the intervention zone in red
  geom_sf(data = intervention_sf, fill = "red", alpha = 0.5, color = "black", size = 1) +
  # Plot the points, color them based on the zone
  geom_sf(data = points_sf, aes(color = zone), size = 1, alpha = 0.7) +
  # set colors for each zone
  scale_color_manual(values = c("Intervention" = "red", "Spillover" = "blue", "Outside" = "grey")) +
  # set plot coordinates
  coord_sf(xlim = c(0, 1000), ylim = c(0, 1000), expand = FALSE) +
  # set plot theme
  theme_minimal() +
  # set plot labels
  labs(title = "Simulated data colored by intervention/spillover/outside", x = "X Coordinate", y = "Y Coordinate") +
  theme(legend.title = element_blank())  # Remove legend title

Compute Response w/ Various Spillover Methods

### Set Constants ###
# Define constants for outcome computation
alpha <- 0.2
beta <- 1

# max distance
max_distance <- sqrt((133-333)^2 + (866-666)^2)
### Compute Response (no spillover) ###
# Calculate distances between each point and the intervention zone
intervention_distance <- points$intervention_distance

# Compute the outcome value: outcome = alpha + beta * x_i
points$outcome_none <- alpha + beta * points$intervention_zone
### Compute Binary Spillover effect ###
# maximum spillover effect
spillover_binary_max <- 0.8

# compute binary spillover outcomes
spillover_binary <- spillover_binary_max * points$spillover_zone
points$outcome_binary <- points$outcome_none + spillover_binary
### Compute Linear Spillover effect ###
# maximum spillover effect
spillover_linear_max <- 1.0

# Compute linear_spillover for points in the spillover zone
linear_spillover <- ifelse(points$spillover_zone == 1, spillover_linear_max*(1-intervention_distance/max_distance), 0) # Linear decay
points$outcome_linear <- points$outcome_none + linear_spillover # add linear spillover adjusted outcome
### Compute Categorical Spillover Effect ###
# Buffer width is based on dividing the distance between spillover and intervention zones by 4
total_buffer_width <- 866 - 666  # The distance from the intervention zone to spillover boundary
buffer_width <- total_buffer_width / 4  # Divide the total width into 4 segments

# Create buffer zones from the intervention zone outward
segment_1 <- st_difference(st_buffer(intervention_sf, buffer_width), intervention_sf)
segment_2 <- st_difference(st_buffer(intervention_sf, buffer_width * 2), st_buffer(intervention_sf, buffer_width))
segment_3 <- st_difference(st_buffer(intervention_sf, buffer_width * 3), st_buffer(intervention_sf, buffer_width * 2))
segment_4 <- st_difference(st_buffer(intervention_sf, total_buffer_width), st_buffer(intervention_sf, buffer_width * 3))


# Loop through each segment and assign the appropriate spillover effect
categorical_spillover <- rep(0, nrow(points))
categorical_spillover[st_intersects(points_sf, segment_4, sparse = FALSE)] <- 0.2
categorical_spillover[st_intersects(points_sf, segment_3, sparse = FALSE)] <- 0.4
categorical_spillover[st_intersects(points_sf, segment_2, sparse = FALSE)] <- 0.6
categorical_spillover[st_intersects(points_sf, segment_1, sparse = FALSE)] <- 0.8

# add effect of 0.25 to any point in spillover zone that hasn't been assigned effect (corners that are left out)
spillover_effect_points <- st_intersects(points_sf, spillover_sf, sparse = FALSE)  # Identify points in spillover zone
no_effect_points <- is.na(categorical_spillover) | categorical_spillover == 0  # Identify points without an effect
categorical_spillover[spillover_effect_points & no_effect_points] <- 0.2

# Now calculate the overall outcome with the spillover effect
points$outcome_categorical <- points$outcome_none + categorical_spillover

# Plot categorical spillover zones (corners included with outermost zone)
# ggplot() +
#   geom_sf(data = spillover_sf, fill = "blue", alpha = 0.2) +  # Spillover zone in light blue
#   geom_sf(data = intervention_sf, fill = "red", alpha = 0.5) +  # Intervention zone in red
#   geom_sf(data = segment_1, fill = "yellow", alpha = 0.3) +  # First buffer zone
#   geom_sf(data = segment_2, fill = "green", alpha = 0.3) +   # Second buffer zone
#   geom_sf(data = segment_3, fill = "orange", alpha = 0.3) +  # Third buffer zone
#   geom_sf(data = segment_4, fill = "purple", alpha = 0.3) +  # Fourth buffer zone
#   theme_minimal() +
#   ggtitle("Concentric Buffer Zones to Spillover Zone Corners") +
#   coord_sf(xlim = c(0, 1000), ylim = c(0, 1000))
### Compute Exponential Spillover effect ###
# constants
spillover_exponential_max <- 1.0 # maximum spillover effect
k <- 5  # Decay rate

# Compute exponential_spillover for points in the spillover zone
exponential_spillover <- ifelse(points$spillover_zone == 1,
                                # exponential decay
                                spillover_exponential_max*exp(-k*(intervention_distance/max_distance)), 0)

# Compute the outcome with the exponential_spillover effect
points$outcome_exponential <- points$outcome_none + exponential_spillover
### Compute Gaussian Spillover effect ###
sigma <- max_distance / 3  # The standard deviation for the Gaussian decay, scaled to the distance

# Compute gaussian_spillover for points in the spillover zone
gaussian_spillover <- ifelse(points$spillover_zone == 1,
                                    # Gaussian decay
                                    exp(-(intervention_distance^2) / (2*sigma^2)), 0)

# Compute the outcome with the gaussian_spillover effect
points$outcome_gaussian <- points$outcome_none + gaussian_spillover
### Compute IDW Spillover effect ###
epsilon <- 0.001  # Small constant to prevent division by zero

# Compute IDW_spillover for points in the spillover zone
idw_spillover <- ifelse(points$spillover_zone == 1,
                        # IDW decay: 1 / (distance_to_intervention + epsilon)
                        1 / (intervention_distance + epsilon), 0)

# Normalize IDW_spillover to ensure it goes to 0 at the spillover boundary
idw_spillover <- ifelse(intervention_distance <= max_distance,
                        idw_spillover * (max_distance - intervention_distance) / max_distance, 0)

# adjust for boundary cases >1.0
idw_spillover <- ifelse(idw_spillover >= 1.0, 1.0, idw_spillover)

# Compute the outcome with the idw_spillover effect
points$outcome_idw <- points$outcome_none + idw_spillover
### Print preview of data ###
head(points)
## # A tibble: 6 × 13
##       x     y intervention_distance intervention_zone[,1] spillover_zone[,1]
##   <int> <int>                 <dbl>                 <dbl>              <dbl>
## 1   578   614                   0                       1                  0
## 2   549   483                   0                       1                  0
## 3   557   451                   0                       1                  0
## 4   700   670                  34.2                     0                  1
## 5   255   992                 335.                      0                  0
## 6   913   613                 247                       0                  0
## # ℹ 8 more variables: zone <chr[,1]>, outcome_none <dbl[,1]>,
## #   outcome_binary <dbl[,1]>, outcome_linear <dbl[,1]>,
## #   outcome_categorical <dbl[,1]>, outcome_exponential <dbl[,1]>,
## #   outcome_gaussian <dbl[,1]>, outcome_idw <dbl[,1]>

Plot response w/ each spillover effect over distance from intervention zone boundary

### Spillover Method Formulas ###
# Define distance range from 0 to 250
distance <- seq(0, 250, length.out = 500)

# Define spillover effects as functions of distance
# binary
binary_spillover <- ifelse(distance >= 0 & distance <= 200, 0.8, 0)

# linear
linear_spillover <- pmax(1 - (distance / 200), 0)

# categorical
categorical_spillover <- ifelse(distance <= 50, 0.8,
                                 ifelse(distance <= 100, 0.6,
                                        ifelse(distance <= 150, 0.4, 
                                               ifelse(distance <= 200, 0.2, 0))))

# exponential
exponential_spillover <- 1.0 * exp(-5 * (distance / 250))

# gaussian
gaussian_spillover <- exp(-(distance^2) / (2 * ((200/3)^2)))

# IDW
idw_spillover <- 1 / (distance + 0.001) * (1 - distance / 250)

### Plot ###
# Create a data frame for plotting
plot_data <- data.frame(
  distance = distance,
  binary_spillover = binary_spillover,
  linear_spillover = linear_spillover,
  categorical_spillover = categorical_spillover,
  exponential_spillover = exponential_spillover,
  gaussian_spillover = gaussian_spillover,
  idw_spillover = idw_spillover
)

# Reshape the data for ggplot
plot_data_long <- melt(plot_data, id.vars = "distance", variable.name = "spillover_method", value.name = "spillover_effect")

# Create the line plot
ggplot(plot_data_long, aes(x = distance, y = spillover_effect, color = spillover_method)) +
  geom_line(size = 1) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.1)) +
  labs(
    title = "Spillover Effect by Distance from Intervention Zone Boundary",
    x = "Distance from Intervention Zone Boundary",
    y = "Spillover Effect Proportion",
    color = "Spillover Method"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

Plot Simulated Data for each of 6 spillover methods

### No Spillover ###
spillover_method <- points$outcome_none

plot_none <- ggplot() +
  # Plot intervention zone boundary
  geom_sf(data = intervention_sf, fill = NA, alpha = 0.5) +
  # Plot spillover zone boundary
  geom_sf(data = spillover_sf, fill = NA, alpha = 0.3) +
  # Plot points
  geom_point(data = points, 
             aes(x = x, y = y, 
                 color = ifelse(intervention_zone == 1, "red", "blue"),  # Color points based on zone
                 alpha = spillover_method),  # Adjust opacity based on response
             show.legend = FALSE) +  # Hide legend for points
  # Add labels
  labs(title = "Simulated Data with No Spillover Effect", x = "X Coordinate", y = "Y Coordinate") +
  # Set theme
  theme_minimal() +
  # Adjust color scale (maintain transparency)
  scale_color_identity() +  # Use specified colors directly
  scale_alpha(range = c(0.1, 0.5))  # Set alpha scale for transparency

# print plot
plot_none

### Binary Spillover ###
spillover_method <- points$outcome_binary

plot_binary <- ggplot() +
  # Plot intervention zone boundary
  geom_sf(data = intervention_sf, fill = NA, alpha = 0.5) +
  # Plot spillover zone boundary
  geom_sf(data = spillover_sf, fill = NA, alpha = 0.3) +
  # Plot points
  geom_point(data = points, 
             aes(x = x, y = y, 
                 color = ifelse(intervention_zone == 1, "red", "blue"),  # Color points based on zones
                 alpha = spillover_method),  # Adjust opacity based on response
             show.legend = FALSE) +  # Hide legend for points
  # Add labels
  labs(title = "Simulated Data with Binary Spillover Effect", x = "X Coordinate", y = "Y Coordinate") +
  # Set theme
  theme_minimal() +
  # Adjust color scale (maintain transparency)
  scale_color_identity() +  # Use specified colors directly
  scale_alpha(range = c(0.1, 0.5))  # Set alpha scale for transparency

# print plot
plot_binary

### Linear Spillover ###
spillover_method <- points$outcome_linear

plot_linear <- ggplot() +
  # Plot intervention zone boundary
  geom_sf(data = intervention_sf, fill = NA, alpha = 0.5) +
  # Plot spillover zone boundary
  geom_sf(data = spillover_sf, fill = NA, alpha = 0.3) +
  # Plot points
  geom_point(data = points, 
             aes(x = x, y = y, 
                 color = ifelse(intervention_zone == 1, "red", "blue"),  # Color points based on zones
                 alpha = spillover_method),  # Adjust opacity based on response
             show.legend = FALSE) +  # Hide legend for points
  # Add labels
  labs(title = "Simulated Data with Linear Spillover Effect", x = "X Coordinate", y = "Y Coordinate") +
  # Set theme
  theme_minimal() +
  # Adjust color scale (maintain transparency)
  scale_color_identity() +  # Use specified colors directly
  scale_alpha(range = c(0.1, 0.5))  # Set alpha scale for transparency

# print plot
plot_linear

### Categorical Spillover ###
spillover_method <- points$outcome_categorical

plot_categorical <- ggplot() +
  # Plot intervention zone boundary
  geom_sf(data = intervention_sf, fill = NA, alpha = 0.5) +
  # Plot spillover zone boundary
  geom_sf(data = spillover_sf, fill = NA, alpha = 0.3) +
  # Plot points
  geom_point(data = points, 
             aes(x = x, y = y, 
                 color = ifelse(intervention_zone == 1, "red", "blue"),  # Color points based on zones
                 alpha = spillover_method),  # Adjust opacity based on response
             show.legend = FALSE) +  # Hide legend for points
  # Add labels
  labs(title = "Simulated Data with Categorical Spillover Effect", x = "X Coordinate", y = "Y Coordinate") +
  # Set theme
  theme_minimal() +
  # Adjust color scale (maintain transparency)
  scale_color_identity() +  # Use specified colors directly
  scale_alpha(range = c(0.1, 0.5))  # Set alpha scale for transparency

# print plot
plot_categorical

### Exponential Spillover ###
spillover_method <- points$outcome_exponential

plot_exponential <- ggplot() +
  # Plot intervention zone boundary
  geom_sf(data = intervention_sf, fill = NA, alpha = 0.5) +
  # Plot spillover zone boundary
  geom_sf(data = spillover_sf, fill = NA, alpha = 0.3) +
  # Plot points
  geom_point(data = points, 
             aes(x = x, y = y, 
                 color = ifelse(intervention_zone == 1, "red", "blue"),  # Color points based on zones
                 alpha = spillover_method),  # Adjust opacity based on response
             show.legend = FALSE) +  # Hide legend for points
  # Add labels
  labs(title = "Simulated Data with Exponential Spillover Effect", x = "X Coordinate", y = "Y Coordinate") +
  # Set theme
  theme_minimal() +
  # Adjust color scale (maintain transparency)
  scale_color_identity() +  # Use specified colors directly
  scale_alpha(range = c(0.1, 0.5))  # Set alpha scale for transparency

# print plot
plot_exponential

### Gaussian Spillover ###
spillover_method <- points$outcome_gaussian

plot_gaussian <- ggplot() +
  # Plot intervention zone boundary
  geom_sf(data = intervention_sf, fill = NA, alpha = 0.5) +
  # Plot spillover zone boundary
  geom_sf(data = spillover_sf, fill = NA, alpha = 0.3) +
  # Plot points
  geom_point(data = points, 
             aes(x = x, y = y, 
                 color = ifelse(intervention_zone == 1, "red", "blue"),  # Color points based on zones
                 alpha = spillover_method),  # Adjust opacity based on response
             show.legend = FALSE) +  # Hide legend for points
  # Add labels
  labs(title = "Simulated Data with Gaussian Spillover Effect", x = "X Coordinate", y = "Y Coordinate") +
  # Set theme
  theme_minimal() +
  # Adjust color scale (maintain transparency)
  scale_color_identity() +  # Use specified colors directly
  scale_alpha(range = c(0.1, 0.5))  # Set alpha scale for transparency

# print plot
plot_gaussian

### IDW Spillover ###
spillover_method <- points$outcome_idw

plot_idw <- ggplot() +
  # Plot intervention zone boundary
  geom_sf(data = intervention_sf, fill = NA, alpha = 0.5) +
  # Plot spillover zone boundary
  geom_sf(data = spillover_sf, fill = NA, alpha = 0.3) +
  # Plot points
  geom_point(data = points, 
             aes(x = x, y = y, 
                 color = ifelse(intervention_zone == 1, "red", "blue"),  # Color points based on zones
                 alpha = spillover_method),  # Adjust opacity based on response
             show.legend = FALSE) +  # Hide legend for points
  # Add labels
  labs(title = "Simulated Data with IDW Spillover Effect", x = "X Coordinate", y = "Y Coordinate") +
  # Set theme
  theme_minimal() +
  # Adjust color scale (maintain transparency)
  scale_color_identity() +  # Use specified colors directly
  scale_alpha(range = c(0.1, 0.5))  # Set alpha scale for transparency

# print plot
plot_idw

Compute Spatial Autocorrelation coefficient